Genome-wide analysis of long noncoding RNAs as cis-acting regulators of transcription factor-encoding genes in IgA nephropathy

Background IgA nephropathy (IgAN) is the most common form of primary glomerulonephritis in the world, but the disease pathogenesis noncoding is yet to be elucidated. Previous studies have revealed regulatory functions for long noncoding RNA (lncRNA) in various diseases; however, the roles of lncRNA in IgAN and regulation of transcription factors (TFs) have been scarcely investigated. Methods Renal tissue samples (n = 5) from patients with IgAN and control samples (n = 4) were collected and RNA sequencing (RNA-seq) was performed. Four software programs were employed for lncRNA prediction. GO (Gene Ontology)/KEGG (Kyoto Encyclopedia of Genes and Genomes) were employed for analysis of the identified differentially expressed genes (DEGs). A regulatory network model of DE lncRNA-TF-DEG was developed, and the levels of expression of key lncRNAs, TFs, and corresponding target genes were assessed using qRT-PCR and immunofluorescence. Results The current study identified 674 upregulated and 1,011 downregulated DE mRNAs and 260 upregulated and 232 downregulated DE lncRNAs in IgAN samples compared with control samples. The upregulated DE mRNAs showed enrichment in cell adhesion and collagen glial fiber organization pathways. The DE lncRNAs-DE mRNAs showing co-expression are associated with transmembrane transport. A novel regulatory network model of lncRNA-TF-DEG has been developed. This study identified seven TFs that are cis-regulated by 6 DE lncRNAs, and show co-expression with 132 DEGs (correlation coefficient ≥ 0.8, P ≤ 0.01), generating 158 pairs that showed co-expression. The lncRNAs NQO1-DT and RP5-1057120.6 were found to be highly expressed in IgAN samples. The TFs vitamin D Receptor (VDR) and NFAT5, along with their target genes were also aberrantly expressed. Conclusion Key lncRNAs and TFs centrally associated with IgAN have been identified in this study. A regulatory network model of lncRNA-TF-mRNA was constructed. Further studies on the genes identified herewith could provide insight into the pathogenesis of IgAN.


Introduction
Immunoglobulin A nephropathy (IgAN) is one of the most common forms of primary glomerulonephritis worldwide and is chiefly characterized by the deposition of IgA in the glomerular mesangium accompanied by increase in mesangial matrix and cellularity and deposition of IgA-immune complexes [1][2][3].Approximately 15%-40% of patients with IgAN typically progress to end-stage renal disease (ESRD) within 20 years of disease onset [4].Although IgAN was first reported over 50 years ago, the pathophysiology of the disease is yet to be elucidated [1,5].
Long noncoding RNAs (lncRNAs) are a type of RNA with length >200 nucleotides that do not typically code for proteins [6].These lncRNAs regulate transcription by altering the chromatin state, regulating transcription factors (TFs), or modulating small host RNAs [7].Several studies have reported the important role of lncRNAs in IgAN.The lncRNA colorectal neoplasia differentially expressed(CRNDE) has been shown to accelerate IgAN development through activation of NLRP3 inflammasome in macrophages [8], whereas lncRNA PTTG3P induces aberrant glycosylated IgA1 production and abnormal B-cell proliferation in IgAN [9].Transcriptional dysregulation is a key mediator of the pathophysiology of various human diseases, such as diabetes, inflammatory disorders, cardiovascular diseases, and several types of cancer [10].Transcriptional dysregulation is associated with intracapillary proliferation in the kidneys of patients with IgAN and has been identified as a risk factor for poor prognosis [11].The transcription factor NF-κB contributes to renal tissue damage by activating the expression of transcriptional regulators [12,13] and is associated with poor prognosis in patients with IgAN [12].The use of GEO database to identify differentially expressed TFs associated with inflammation could reveal potential targets for the treatment of IgAN [14].Therefore, both lncRNAs and TFs play important roles in the pathogenesis of IgAN.A few studies have illustrated the role of lncRNA in the development and progression of other diseases through regulation of TFs.Cooperative regulation among two lncRNAs (MALAT1 and NEAT1), three TFs (NF-κB, NFE2L2, and PPARG), and a set of miRNAs has been shown to regulate downstream differentially expressed genes (DEGs) [15].The lncRNA MAGI2-AS3 has been shown to regulate ACY1 through interaction with the transcription factor HEY1 to prevent tumor growth and angiogenesis in clear cell renal cell carcinoma [16].To date, such a regulatory association between lncRNA and transcription factors (TFs) has not been demonstrated in IgAN expression.We hypothesize that differential lncRNA expression of IgAN could regulate differential expression of TFs and thereby, differential expression of downstream genes.The present study aims to ascertain a role for lncRNAs in transcriptional regulation of IgAN through the identification of genome-wide lncRNA-TF-DEG network.

Tissue specimens
Renal tissue samples from patients with primary IgAN (IgAN; n = 5) and histologically normal kidney tissue 5 cm adjacent to the renal cancer tissue (control; n = 4) were collected between February 1, 2021, and May 1, 2021, in the First Affiliated Hospital of Zhengzhou University (five IgAN and four control samples were subjected to RNA-seq, whereas samples from another three IgAN patients and another three control samples were employed for validation).The inclusion criteria for IgAN were as follows: 1) primary IgAN diagnosed via kidney biopsy in the Renal Department of the First Affiliated Hospital, Zhengzhou University, 2) complete clinical and pathological record/data for the recruited patients, 3) glomerular number >10 in the renal tissue sample, and 4) glucocorticoids or other immunosuppressors not administered to the patients.All tissue samples were immediately frozen in liquid nitrogen and stored a t −80˚C for RNA extraction.The study protocol was approved by the Institutional Ethics Board of the First Affiliated Hospital of Zhengzhou University and the written informed consent was obtained from all the participants.

Reads alignment and DEG analysis
Clean reads were aligned with the human GRCh38 genome using HISAT2 [17].Uniquely mapped reads were used to calculate read numbers and fragments per kilobase of exon per million mapped fragments (FPKM) for each gene.The expression levels of genes were evaluated using FPKM.Analysis of differential expression of genes was carried out using the software DEseq2 and determined via fold change (FC) and false discovery rate (FDR), with FC � 2 or � 0.5 and FDR < 0.05 indicating significant differential expression.

Prediction of lncRNA
Four software programs, namely, CPC2 [18], LGC [19], CNCI [20], and CPAT [21].were employed for lncRNA prediction.Of the noncoding transcripts identified using the above programs, the following transcripts were discarded: those showing overlap with known coding sequences, those <200 nucleotides in length, those having a potential coding ability, or those encoded by a sequence located at a distance of <1,000 bp from the nearest gene as observed from assembly results.The predicted lncRNA obtained using a Venn-diagram intersection of lncRNA transcripts generated using the four software programs mentioned above were employed for subsequent analysis and processing.

Identification of cis-regulatory targets of lncRNA
The threshold of co-location was set as 100 kb both upstream and downstream of lncRNA in the trans-regulatory relationship pair [22], followed by screening of the lncRNA-target mRNA relationship pairs for co-expression using the criteria of absolute correlation values >0.6 and P � 0.01.The intersection of the two data sets of co-location and co-expression was employed to obtain the cis targets of lncRNA.The lncRNA cis-targeted TFs (cis-TFs) were subsequently filtered out from the above by referring to a list of 1,796 TFs (HumanTFDB; http://bioinfo.life.hust.edu.cn/HumanTFDB#!/).

Functional enrichment analysis
Gene Ontology (GO) terms and KEGG pathways were identified using KOBAS 2.0 server to sort out functional categories of DEGs [23].Hypergeometric test and Benjamini-Hochberg procedure for FDR control were used to define enrichment of each term.Reactome (http:// reactome.org)pathway profiling was also employed for functional enrichment analysis of the sets of selected genes.

Quantitative real-time PCR (qRT-PCR) analysis
Total RNA was extracted from renal tissue samples using TRIzol reagent (CWBIO, CHINA) and used as template for cDNA synthesis using cDNA Synthesis Kit (US EVERBRIGHT, CHINA).Subsequently, qRT-PCR analysis was performed using SYBR-Green Master Mix (US EVERBRIGHT, CHINA) and AI Q3 machine (Thermo Fisher Scientific, USA).Relative mRNA levels in IgAN vs control tissue samples normalized to GAPDH were calculated, and the relative change in gene expression was analyzed using the 2−ΔΔCT method.The primers used are shown in Table 1.

Immunofluorescence
Tissue sections were prepared as detailed previously using fixing, dehydration, slicing, baking, dewaxing, and hydration techniques.The sections were repaired by heating and boiling for 20 min in EDTA-antigen retrieval buffer (pH 8.0).After cooling, the sections were washed thrice with PBS solution (pH 7.4) in a rocker for 5 min each, followed by the dropwise addition of blocking solution (3% BSA) and incubation at 37˚C for 30 min.The primary antibodies rabbit anti-NFAT5 (ab3446, Abcam, England) and rabbit anti-vitamin D Receptor (VDR) (GB114918, Servicebio, China) were added dropwise and the sections were incubated overnight at 4˚C in a wet box.The sections were subsequently washed three times with PBS (pH 7.4) as described, and incubated with biotinylated goat anti-rabbit Immunoglobulin G (IgG) for 50 min at room temperature under dark conditions.This was followed by washing three times with PBS (pH 7.4) as above, and subsequent incubation with DAPI solution at room temperature for 10 min under dark conditions.The sections were again washed three times with PBS as above, followed by incubation with spontaneous fluorescence quenching reagent for 5 min.This was followed by washing under running water for 10 min, and use of anti-fade mounting medium for mounting.The sections were observed under a fluorescence microscope and images were obtained.

Western blot
Total protein of tissue were extracted using a total protein extraction buffer (Beyotime, China).We performed western blotting according to the traditional western blot assay method.β-ACTIN was utilized as an internal control, Antibodies (1:1000) against VDR was purchased from Servicebio (Wuhan, China).

Statistical analysis
Principal component analysis (PCA) was performed using R package factoextra (https://cloud.r-project.org/package=factoextra)to show the clustering of samples in different groups.After normalizing the reads by Tags Per Million for each gene, an in-house script (sogen) was used for visualizing next-generation sequence data and for genomic annotations.The heatmap package (https://cran.rproject.org/web/packages/pheatmap/index.html) in R was used to perform clustering based on Euclidean distance.Student's t-test was used to conduct comparisons between the two groups, considering a p-value of < 0.05 as statistically significant.

Identification of lncRNAs expressed in IgAN
The design and workflow of the present study, which utilized a total of 9 renal tissue samples (IgAN, n = 5; control, n = 4) for whole-genome RNA-seq analysis, are presented in Fig 1A .Of the previously identified lncRNA genes, a total of 9,884 lncRNAs were observed in IgAN samples as opposed to 9,442 lncRNAs in the control, including 8,517 lncRNAs with overlapping expression.In addition, 578 newly predicted lncRNAs in IgAN samples and 549 lncRNAs in control samples were identified using four methods (CPC2, LGC, CNCI, and CPAT) and predicting the overlapped LncRNAs (S1 Fig in S1 File), 457 overlapped LncRNAs were included (Fig 1B and 1C).Furthermore, PCA for all mRNA-and lncRNA-encoding genes in IgAN and control samples revealed significant differences between the two groups (Fig 1D and 1E).

Analysis of differential expression of lncRNAs and mRNAs in IgAN samples
Given the differences between IgAN and control renal tissue, PCA was employed to analyze the mRNAs (DE mRNA) and lncRNAs (DE lncRNA) that were differentially expressed between the two samples.

TFs under cis-regulation by DE lncRNAs in IgAN samples
The scatter plot in Fig 3A shows the log-transformed FC of DE lncRNA and cis-regulated genes between IgAN and control samples.Cis-regulation is an important mechanism of regulation of gene expression by lncRNAs.In order to ascertain the regulatory functions of DE lncRNAs on genes 10 kb downstream, co-expression analysis was carried out on DE lncRNAs

Differential expression of lncRNA, cis-TF, and target DEGs in IgAN and control samples
The

Validation of lncRNAs and gene regulatory network
The expression levels of lncRNAs NQO1-DT and RP5-1057120.6 and the downstream TFs and targets were validated using qRT-PCR in IgAN and control samples from another cohort (Fig 5).Immunofluorescence analysis was also employed to validate the expression levels of the TFs (Fig 6).As expected, the expression levels of lncRNA NQO1-DT and lncRNA RP5-1057120.6 were significantly higher in the IgAN compared with the control samples.Conversely, IgAN samples exhibited much higher levels of NFAT5 but lower levels of VDR compared with control samples.The differences in expression levels of PAX8, PCDH7, ADM, and SLIT2 between IgAN and control samples were consistent with previous results; statistically significant differences (P>0.05) were not observed for PODXL and SLC22A8, whereas WT1 and GAS1 demonstrated an opposite trend compared with the above results.

Discussion
IgAN is one of the most common forms of primary glomerulonephritis globally, with approximately 15%-40% of patients progressing to ESRD within 10-20 years of disease onset [4].However, the current understanding of the disease pathophysiology remains incomplete.The mammalian genome encodes numerous lncRNAs, which play key functional roles in various biological processes [25,26].In the present study, we identified DE lncRNAs (260 upregulated and 232 downregulated) and DE mRNAs (674 upregulated and 1,011 downregulated) that show differential expression in IgAN compared with control samples.The heatmap of 30 DE lncRNAs that demonstrated highest upregulation or downregulation in IgAN samples revealed that the downregulated lncRNA RPPH1 and upregulated lncRNA CRNDE are particularly significant.Shen et al. previously reported that the lncRNA CRNDE accelerates IgAN progression by promoting NLRP3 inflammasome activation in macrophages [8], whereas Zhang et al. reported that lncRNA RPPH1 promotes inflammatory response and cell proliferation in mesangial cells in diabetic nephropathy [27], but its role in IgA nephropathy is unclear.The GO and KEGG pathway analyses carried out in the present study revealed that the upregulated DE mRNAs were mainly enriched in cell adhesion and collagen glial fiber organization pathways.Previous studies have also demonstrated roles for cell adhesion and collagen fibril pathways in IgAN [28][29][30][31].Cell adhesion and collagen fibrils are central to the infiltration of inflammatory cells into the renal interstitium, leading to tubulointerstitial inflammation and subsequent interstitial fibrosis in various kidney diseases.These observations suggest that lncRNAs could play a crucial role in IgAN and may correlate with disease pathogenesis.
The present study has established a regulatory network consisting of lncRNA and TFs-DEG, whereby lncRNAs regulate TFs, which subsequently yield changes in expression of target genes.A total of 7 TFs under cis-regulation of 6 DE lncRNAs were identified in the present study, and these were observed to be co-expressed with 132 DEGs (correlation coefficient � 0.8, P � 0.01), accounting for a total of 158 pairs that showed co-expression.GO analysis revealed enrichment of these 7 cis-TFs in transcriptional regulation pathways, and of the DEGs targeted by these cis-TFs in pathways regulating cell proliferation and apoptosis.Zhang et al. reported that the deposition of IgA1-IgG immune complexes in the mesangial region of the glomeruli in IgAN stimulates mesangial cell proliferation and increases the secretion of cytokines, chemokines, and extracellular matrix proteins; moreover, podocyte apoptosis has been observed in a mouse model of IgAN, which is consistent with our results [32].
Previous studies have revealed important roles for two TFs, VDR and NFAT5 [24], and their downstream target genes in the cis lncRNA-mRNA network.The transcription factors VDR and Nuclear Factor of Activated T Cells 5 (NFAT5) along with their corresponding downstream targets, were found to play an important role in the cis lncRNA-TF-DEG regulatory network.Adrenomedullin (ADM), Wilms' tumor suppressor 1 (WT1), Slit guidance ligand 2 (SLIT2), Growth arrest-specific protein-1 (GAS1), Solute carrier family 22 member 8 (SLC22A8), and Podocalyxin-like protein (PODXL) are the target genes of the transcription factor VDR and exhibit co-expression with it, whereas Paired box gene 8 (PAX) and Protocadherin 7 (PCDH7) are the co-expressed target genes of NFAT5.Notable differences were observed in the expression of the TFs and target genes between IgAN and control samples.These differences were validated using qRT-PCR; statistically significant differences were observed in the expression of TFs VDR and NFAT5 and the corresponding downstream target genes between IgAN and control samples.This in turn indicates a key role for the TFs VDR and NFAT5 in the pathogenesis of IgAN.But in the validation experiments, we found the genes are not consistently changed.We speculated that there may be several reasons for this result.(1) We chose different tissues for validation, there may be individual differences among different individuals.(2) Heterogeneity of IgA nephropathy.
Mo et al. reported higher levels of blood urea nitrogen, serum phosphorus, mesangial cellularity/cell proliferation, interstitial fibrosis/tubular atrophy and crescents in VDR rs2228570T allele carriers compared with non-T allele carriers, whereas eGFR and 25-Hydroxyvitamin D3 were lower in T-allele carriers compared with non-T allele carriers among IgAN patients, indicating a role for VDR in the pathogenesis of IgAN [33].ADM, a vasodilator peptide abundantly expressed in kidneys, shows enrichment with pathways that positively regulate apoptosis.Previous studies have reported that decrease in glomerular ADM in IgAN patients occurred prior to the onset of renal dysfunction and could be associated with the proliferation of glomerular mesangial cells [34].SLIT2 showed enrichment with pathways involved in negative regulation of cell proliferation and positive regulation of apoptosis.The SLIT2/ROBO2 signaling pathway was found to inhibit non-muscle myosin IIA activity and destabilize kidney podocyte adhesion [35].
PAX8 and PCDH7 are both target genes of NFAT5 that show co-expression with the TF.PAX8 has been reported to show enrichment with pathways that regulate apoptosis and RNA polymerase II transcription, and the PAX8 protein is believed to play a role in regeneration and recovery following acute kidney injury [36].PCDH7 is known to play a role in platelet degranulation and cell adhesion.Sethi et al. reported that 5.7% of patients with PLA2R-negative membranous nephropathy were positive for PCDH7 staining within glomeruli, indicating that it could be a potential biomarker for a distinct, previously unidentified type of membranous nephropathy [37].
Of course, there are still several limitation in our study.First, we only verified that lncRNA, transcription factor and mRNA were differentially expressed between IgA nephropathy and normal control, but did not explore the specific mechanism.Second, the results need to be verified in more samples.The DEGs and transcription factors screened in this study were mainly expressed in the renal tubulointerstitial section, so in the subsequent studies, we can mainly focus on the function of the target genes on renal tubular and interstital.

Conclusion
A novel function for dysregulated lncRNA in regulating the differential expression of TFs and thereby altering expression of downstream genes was identified in the present study.A regulatory network model of cis DElncRNAs-TF-DEG was developed, which may promote understanding of the pathogenesis of IgAN.

Fig 1 .
Fig 1. Identification of expressed lncRNAs in IgA nephropathy.A. Workflow of lncRNA prediction and analysis in this research.B-C.Venn diagram showing the number of detected lncRNAs.The known and new predicted lncRNAs were detected.(Expressed at least 1 sample with FPKM > 0) D. Principal component analysis (PCA) based on FPKM value of all mRNA.The ellipse for each group is the confidence ellipse.E. Principal component analysis (PCA) based on FPKM value of all lncRNAs.The ellipse for each group is the confidence ellipse.https://doi.org/10.1371/journal.pone.0304301.g001

Fig 5 .
Fig 5. Relative expression levels of lncRNAs, transcription factors (TFs) and mRNAs in validation cohort.(The experiments were repeated three times.Student's t test was used to analyze the statistical significances of differences between two groups).https://doi.org/10.1371/journal.pone.0304301.g005

Fig 6 .
Fig 6.Transcription factors VDR and NFAT5 is different expression in IgAN and healthy control.A. Immunofluorescence analysis showing that VDR expression level in healthy control (HC) was higher than in IgAN.B. Immunofluorescence analysis showing that NFAT5 expression level in IgAN was higher than in healthy control.https://doi.org/10.1371/journal.pone.0304301.g006